function [statedr,etadr,stt_o,eta_o,smooth_o]=kalman_dr1_loop(y,xdata,Z,Cobs,T,R,Qmat,A,addin,...
    pos_stdr,flag_robust)
% ==================================================================
% [stt_o,eta_o,smooth_o,ptt_o,state_dr,eta_dr,yfor_o,yferr_o]=kalman_dr1(y,
% xdata,Z,Cobs,T,R,Qmat,A,addin)
%
% KALMAN_DR1_LOOP
% 
% This is the same as KALMAN_DRALL but for a single draw 
% Reported state draw is trimmed to just the column positions passed in 
% POS_STDR. If that argument is missing or [], then all states are
% returned 
% Use this code with LOOP over DRAWS when only a few states are of
% interest 
%
% See also KALMAN_DR1 which does not require the indicator and uses a
% different subfunction
%
% Checked against KALMAN_DRALL.m in EP_REALRATE. 
% Alejandro Justiniano June 26 2009
% =============================================================
if nargin < 10 || isempty(pos_stdr)
    pos_stdr=(1:size(T,1)); 
end 

if nargin < 11; flag_robust=0;end; 
if flag_robust~=0 && flag_robust~=1;error('Flag robust must be 0 or 1');end 

% 1. Initial check for dimension
% -------------------------------
[nobs,nz]=size(y);
[ny,nx] = size(R);
ns=size(T,1);
Q=Qmat'*Qmat;
RQ=R*(Qmat'); 

% 2. Removal of Auxiliary data and means
% ydem
% ---------------------------------------
if ~isempty(xdata);
    %flag_xdata=1;
    if size(xdata,1)~=nobs;error('X and Y must have same number of rows');end
    if size(xdata,2)~=size(A,2);error('Cols A must match Cols of X');end;
    if size(A,1)~=size(T,1);error('Rows A must match dimension state space');end;
    ydem=y-xdata*(A');
else
    %flag_xdata=0;
    ydem=y;
end

if any(Cobs~=0)==1;
    ydem=ydem-repmat((Z*Cobs)',[nobs 1]);
end;

% 3. Initialization
% ==================
% I. Non-diffuse
% PSHAT P0|0
% SHAT  s0|0 
if addin.flag_pzero==1
    pzero=addin.pzero;
else
    pzero=lyapunov_symm(T,R*Qmat*(Qmat')*(R'));
    pzero=0.5*(pzero+pzero'); 
    if ~isempty(find(isnan(pzero))~=0)
        error('PSHAT contains Nan entries')
    end
end
if addin.flag_szero==1;
    szero=addin.szero;
    if length(szero)~=ns;error('SHAT must be NSx1');end
else
    szero=zeros(ns,1);
end

% ====================================================
% 4. Filter and Smooth YDEM 
% Resulting matrices have _o for orignal 
% =====================================================
Ztr=Z'; %clear Z;
Rtr=R'; %clear R;
[stt_o,eta_o,smooth_o]=kf_forbackfast_sub(ydem',szero,pzero); 
stt_o=stt_o'; 


% ============================================= 
% 6.    Simulate the innovations and the data
% w+
% y+
% s+    alpha+
% =============================================
wplus=randn(nobs,nx)*Qmat;
[yplus,stplus]=simssmodel(wplus,T,R,Z,szero,pzero);
[stt_plus,eta_plus,smooth_plus]=kf_forbackfast_sub(yplus',szero,pzero);
statedr=smooth_o(:,pos_stdr)-smooth_plus(:,pos_stdr)+stplus(:,pos_stdr);
etadr=(eta_o-eta_plus+wplus);

% =========================================================================
% =========================================================================
% =========================================================================
    
    % =====================================================================
    % Soubroutine 
    % KF_FORBACKFAST_SUB 
    % Same as KF_DK but in subroutine mode 
    % STT and PTT DO NOT have the inital S0 and P0 in the first entry 
    % =====================================================================
    function [stt,etamat,smooth_st]=kf_forbackfast_sub(yy,a_zero,p_zero)
    %------------------------------------
    % a.FORWARD FILTER
    % -----------------------------------
    vstar=zeros(nz,nobs);
    finmat=zeros(nz,nz,nobs);
    kpartg=zeros(ns,nz,nobs);
    logLnc=zeros(nobs,1);
    %yfor=zeros(nz,nobs);
    % Matrices with one additional entry (initialization)
    stt=zeros(ns,nobs+1);
    ptt=zeros(ns,ns,nobs+1);
    stt(:,1)  =a_zero;
    ptt(:,:,1)=p_zero;
   
    ii=1; 
    for ii=1:nobs;
        [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vstar(:,ii),finmat(:,:,ii),kpartg(:,:,ii)]...
            =kf_dk(yy(:,ii),Z,stt(:,ii),ptt(:,:,ii),T,RQ);
    end
    stt=stt(:,2:end);
    ptt=ptt(:,:,2:end);
    
    %------------------------------------
    % b.BACKWARD FILTER
    % -----------------------------------
    etamat   =zeros(nx,nobs);
    smooth_st=zeros(ns,nobs);
    
    rstar=zeros(ns,1);
    rmat=zeros(ny,nobs);
    
    [rstar,etamat(:,end)]=smoothdis(rstar,Q,Rtr,Ztr,finmat(:,:,end),zeros(ns),vstar(:,end));
    smooth_st(:,nobs)=stt(:,end)+ptt(:,:,end)*rstar;
    rmat(:,nobs)=rstar;
    
    ii=nobs; 
    for ii=(nobs-1):-1:1;
        [rstar,etamat(:,ii)]=smoothdis(rstar,Q,Rtr,...
            Ztr,finmat(:,:,ii),((T-T*kpartg(:,:,ii)*Z)'),vstar(:,ii));
        smooth_st(:,ii)=stt(:,ii)+ptt(:,:,ii)*rstar;
        rmat(:,ii)=rstar;
    end
    
    % =============================================================
    % c. Check if flag_robust==1
    % Perform another disturbance smoother and check can
    % recover the observables
    % ==============================================================
    if flag_robust==1
        smooth_st2=zeros(ns,nobs);
        smooth_st2(:,1)=T*a_zero+p_zero*rstar;
        ycheck=zeros(size(yy));
        ycheck(:,1)=Z*smooth_st2(:,1);
        for ii=2:nobs;
            smooth_st2(:,ii)=T*smooth_st2(:,ii-1) + R*(etamat(:,ii));
            ycheck(:,ii)=Z*smooth_st2(:,ii);
        end
        ydiscrep=max((abs(ycheck'-yy')));
        if max(ydiscrep > 1e-4);
            error('SMOOTH innovations do not recover observables');
        end
        
        % =============================================
        % Check also difference between smooth states
        if max(max(abs( smooth_st2-smooth_st ))) > 1e-4
            error('SMOOTH states do not match')
        end
    end
    
    etamat=etamat';
    smooth_st=smooth_st';
    
    end

end 